using Pkg Pkg.activate(".") using DataFrames, CSV, StatsPlots, Downloads, DataFramesMeta, Dates, Statistics using Turing, Memoization, ReverseDiff, Serialization using LinearAlgebra, Distributions
In late June 2022 the Supreme Court of the United States (SCOTUS) issued a ruling in New York State Rifle and Pistol Association v Bruen which ended the possibility for any state to issue concealed weapons (CCW) permits on the basis of judgement by government employees (known as "may issue" laws). For example in New York an applicant had to convince a police department that they had a "special need" for protection (such as for example a credible death threat) and in California the Sherriff of each county had a similar discretionary ability to decide whether there was a "sufficient cause" for issuing a permit. Simple self defense was not enough.
In the Bruen ruling, either permitless systems as now exist in 25 states, or CCW permit systems in which "objective" standards apply were the only allowable options (known as "shall issue" laws). Objective standards are for example "has never been convicted of a felony" or "has passed a required safety class" and other similar requirements which are factually verifiable by anyone.
After this ruling, as with most rulings liberalizing gun laws, political opponents of the ruling predicted "wild west" style shootouts in the streets. But in fact at the time of the ruling 43 states had either permitless carry (25 states) or a "shall issue" CCW permit system (18 states).
The usconcealedcarry.com website shows the following "shall issue" map. The only misleading data point is that Vermont is shown as red, because they have never regulated concealed carry in any way, it was explicitly constitutionally guaranteed from the beginning and so there is no CCW permit to issue at all.

Let's begin by loading several datasets. We have data on background checks each month, data matching FIPS codes for states, and CDC WONDER data on firearms homicides, as well as all homicides. We separate homicides from firearms suicides as they have a different cause and require different prevention techniques than violent crime. In general firearms suicides are about 2x as big as all homicides.
Loading and munging the data into useful form:
if !Base.Filesystem.ispath("data/nics-firearm-background-checks.csv") Downloads.download("https://github.com/BuzzFeedNews/nics-firearm-background-checks/raw/master/data/nics-firearm-background-checks.csv","./data/nics-firearm-background-checks.csv") end if ~Base.Filesystem.ispath("data/fipscodes.csv") Downloads.download("https://www2.census.gov/geo/docs/reference/state.txt","./data/fipscodes.csv") end fipscodes = @chain CSV.read("data/fipscodes.csv",DataFrame) begin @subset(:STATE .< 60) end # manually downloaded CDC wonder gun homicide data using online form gunhomicidesnew = CSV.read("./data/wonder-gun-homicide-byyear.csv",DataFrame) gunhomicidesold = CSV.read("./data/cdc-wonder-firearm-homicide-1979-1998.csv",DataFrame) gunhomicides = [gunhomicidesnew ; gunhomicidesold] gunhomicides = @chain rename(gunhomicides,Dict("Crude Rate" => "crudedeathrate", "State Code" => "StateCode", "Year" => "year", "Year Code" => "YearCode")) begin @subset(.! ismissing.(:Deaths)) @transform(:crudedeathrate = map(x -> isnothing(x) ? missing : x, tryparse.(Float64,String.(:crudedeathrate))), :state = String31.(:State)) @subset(.! ismissing.(:crudedeathrate)) @orderby(:StateCode,:year) end gunchecks = CSV.read("./data/nics-firearm-background-checks.csv",DataFrame) startpops = @chain gunhomicides, @subset(:year .== 1999), @select(:State, :StartingPop = :Population,:StartDeathRate = :crudedeathrate) gunsalesest = @chain gunchecks, @orderby(:state,:month), @select(:state,:month,:salesest = 1.1 .* (:handgun .+ :long_gun) .+ 2.0 .* :multiple), groupby(:state), @transform(:cumsales = cumsum(:salesest),:year=year.(:month)), @subset(month.(:month) .== 6) joined = leftjoin(gunhomicides,gunsalesest, on=[:State => :state,:year => :year]; matchmissing=:notequal,makeunique=true) joined = @chain leftjoin(joined,startpops; on=[:state => :State],makeunique=true) begin @subset(.! ismissing.(:Deaths) .&& .! ismissing.(:StartingPop)) @subset(:StateCode .< 60) @transform(:cumgunrate = (:cumsales .+ 1.0 .* :StartingPop) ./ :Population ) end joined = @orderby(leftjoin(joined,fipscodes,on = :state => :STATE_NAME),:StateCode,:year)
Although we don't have a proper starting value for the beginning of the period, we know that for quite a while now the number of guns per capita in the US is on the order of 1. We therefore put a starting value of 1 for every state even though in fact this should vary somewhat likely between 0.7 and 1.5 or so, and then cumulatively total up the approximate number of firearms purchased since then and divide by the running population estimate, to calculate a firearms per capita time series.
We then plot firearms homicides vs firearms per capita as a parametric curve with time not visible. If purchasing guns causes gun homicides, then as guns per capita increases through time so will firearms homicides. However there are also other models that could cause such trends.
Nevertheless, if gun purchases do cause gun homicides, we should see these positive trending lines.
display(plot(joined.cumgunrate,joined.crudedeathrate; group=joined.state, legend=false, xlab="Gun Rate (guns/capita)", ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Gun Prevalence", alpha=.5)) display(plot(joined.cumgunrate,joined.crudedeathrate; group=joined.state, legend=false, xlab="Gun Rate (guns/capita)", ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Gun Prevalence", alpha=.5, ylim=(0,12)))
In addition we plot how Homicides have changed through time, and how firearms per capita has changed through time
display(plot(joined.year,joined.crudedeathrate; group=joined.state, legend=false, xlab="Year", ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Time", alpha=.5)) display(plot(joined.year, joined.cumgunrate; group=joined.state, legend=false, xlab="Year", ylab="Firearms Per Capita Estimate (guns/person)", title="Gun Ownership vs Time", alpha=.5, ylim=(0,4),xlim=(1999,2020)))
The most recent trend in firearms carry in the US is the conversion of "shall issue" states to a permitless carry system. In this system everyone over a certain age (18 or 21 depending on state) who is legally able to possess a pistol may carry that pistol concealed in public subject to relatively few restrictions (such as not allowed in secure areas with metal detectors etc). For the moment we ignore the issue of "open carry" which has different issues such as intimidation purposes.
Of all the systems for regulating firearms in public, this is the least restrictive. There are various models one might have of how this would affect gun homicides. Some might argue armed citizens deter crime, others might argue it makes it easier to commit crime, others might argue the effect if any would be minimal as criminals were already carrying guns and only relatively "model citizens" would start carrying, and they are generally not involved in crime in terms of committing them or as victims. It's worthwhile to look at the data to see what hypotheses are comopatible with the observed trends.
We plot firearms homicide rates through time for each state. We place a vertical red line at the point in time where the constitutional carry law was signed if any.
concarry = leftjoin(fipscodes,CSV.read("./data/ConstCarryDates.csv",DataFrame),on = :STUSAB => :StateCode) concarry.date = map(x -> ismissing(x) ? missing : Date(div(x,10000),div(mod(x,10000),100),1),concarry.LawDate) concarry.year = map(x -> ismissing(x) ? missing : year(x),concarry.date) concarry = @subset(concarry,:STATE .<= 59) #ignore minor territories let p = [] for (st,stusab,statename,statens,lawdate,stdate,year) in eachrow(concarry) d = year sub = @subset(joined,:state .== statename) if(nrow(sub) > 0) pp = plot(sub.year,sub.crudedeathrate,ylim=(0,15),size=(500,500), title="$statename\nhomicides",xlab="Year",ylab="Gun Homicide/100k/yr",legend=false) if(!ismissing(d)) pp = plot!([(d,0),(d,15)]) end push!(p,pp) end end display(plot(p...; size = (2000,2000))) end
Looking at the raw data we can see that there are a number of states that have an "accelerating" trend of increased firearms homicide violence in the period from about mid 2014 to 2015 onward. Some of these have converted to Constitutional Carry near the beginning of that period or during that period, with the trend beginning before the carry provision for several of the states. Also many states have passed their laws in the last few years and have no data post-law.
Let's look at these states grouped into geographic regions, because trends in regions may be more obvious. We will start with the following groupings:
stategroups = [["WA","OR","CA","NV","AZ"],["ID","MT","WY","UT","CO","NM","TX"],["ND","SD","NE","KS","OK","IA","MO"], ["MN","WI","IL","IN","MI","OH"],["AR","LA","MS","AL","GA","FL","SC"],["TN","NC","KY","WV","VA","MD","DE"], ["PA","NJ","NY","CT","RI","MA","VT","NH","ME"],["AK","HI","DC"]] stategroupdf = DataFrame(STUSAB = reduce(vcat,stategroups),stgroup = reduce(vcat,[[k for j in 1:length(stategroups[k])] for k in 1:length(stategroups)])) stategroupdf = leftjoin(fipscodes,stategroupdf; on = :STUSAB) @assert(length(unique(reduce(vcat,stategroups))) == 51) # 50 states plus DC let p1 = [], p2=[] for statelist in stategroups sub = joined[in.(joined.STUSAB, Ref(statelist)),:] ccdates = @chain @subset(concarry, in.(concarry.STUSAB,Ref(statelist))) begin @subset(.! ismissing.(:year)) end ccdatesenact = leftjoin(ccdates,sub, on = [:STUSAB => :STUSAB, :year => :year], makeunique = true) ccdatesenact = @subset(ccdatesenact, .! ismissing.(:crudedeathrate)) #display(ccdatesenact) #display(sub) p = plot(sub.year,sub.crudedeathrate,group=sub.STUSAB,ylim=(0,20),legend=:topleft,linewidth = 3, thickness_scaling=1, ylab="Homicides/100k population/yr",xlab="Year") if nrow(ccdatesenact) > 0 p = scatter!(ccdatesenact.year,ccdatesenact.crudedeathrate,markersize=6,group = ccdatesenact.STUSAB) end push!(p1,p) p = plot(sub.cumgunrate,sub.crudedeathrate,markersize=6,linewidth=3,thickness_scaling=1, xlab="Gun Ownership (guns/person)", ylab="Homicide rate (per 100k)", group=sub.STUSAB) push!(p2,p) end display(plot(p1...; size=(1000,1000))) display(plot(p2...; size=(1000,1000))) end
Each state has its own overall level of firearm violence, but often trends in the same direction relative to its region even if its overall level is higher or lower. This suggests a model in which we measure each state using a dimensionless number relative to some baseline level. For example the average rate in the years 1999,2000,2001
baselines = @chain joined begin @subset(in.(:year,Ref((1999,2000,2001)))) @by(:STUSAB, :baseline = mean(:crudedeathrate)) end bljoined = @chain leftjoin(stategroupdf,leftjoin(baselines,joined,on = :STUSAB, makeunique=true),on=:STUSAB,makeunique=true) begin @subset(.! ismissing.(:crudedeathrate)) @orderby(:stgroup,:STUSAB,:year) end function getccd(gr) ccd = innerjoin(gr,@subset(concarry,.! ismissing.(:year)), on = [:STUSAB => :STUSAB, :year => :year], makeunique = true) ccd = @subset(ccd, .! ismissing.(:crudedeathrate)) ccd end function makerelplot(gr) ccd = getccd(gr) # display(ccd) p = plot(gr.year,gr.crudedeathrate ./ gr.baseline,group=gr.STUSAB,legend=:topleft,ylim=(0.0,2.5),xlab="Year",ylab="Relative Rate") if nrow(ccd) == 1 scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline) elseif nrow(ccd) > 1 scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline,group = ccd.STUSAB) end p end
makerelplot (generic function with 1 method)
Each state plotted in a group with other geographically similar states. Firearm homicide rate relative to average recorded rate 1999-2001.
plots = [makerelplot(gr) for gr in groupby(bljoined,:stgroup)] plot(plots... ; size=(1000,1000),linewidth=2,xlab="Year",ylab="Relative Rate")
We can try to estimate the effects of these laws by building a model of the overall process. Because we can't resurrect people from the dead, the homicide rate can never be negative. It makes sense then to model the homicide rate on a logarithmic scale.
We estimate the behavior of the states as an overall level, plus a shape that is common to the region. Individual states are then allowed a small perturbation to the regional shape.
This small perturbation method was the only way that seemed to make sense to estimate the counterfactual. Clearly none of the states will perform exactly like the average. So if we use the regional average, we will bias towards a lot of noise. Obviously, there is no effect if no law was passed, so for those states with no law passed, the counterfactual is the same as the actual, and we then should allow a size of perturbation which allows the function to fit the actuals for those states. By modifying the parameter statecoefpert we can make it as small as possible while still causing all the non-law states to fit reasonably well. At this point whatever lack of fit there is post-law in those states, could be at least partially attributed to the law. The choice of value of log(1.25)/2.0 ~ 0.112 seemed to fit the bill, leaving only Maryland the only non-law not quite well fit by the model. Any smaller value left both MD and NY fitting poorly.
The effect of constitutional carry law would then be estimated as the log firearm homicide rate which exceeds the state level counterfactual estimate in the years after the passage of the law, with a transient onset curve included in the model.
We impose an informal smoothness requirement on the counterfactual estimates by using a compact radial basis function expansion with one center every 4 years, and a maximum radius of 7 years so that there is overlap between adjacent centers.
function bumpfun(x,c,scale) stdx = (x-c)/scale if stdx < -1.0 || stdx > 1.0 0.0 else exp(1.0-1.0/(1.0-stdx^2)) # goes to zero at -1 and 1, and 1 at x=0 end end function timeser(x,coefs,centers,scale) f = 0.0 for (a,c) in Iterators.zip(coefs,centers) f += a*bumpfun(x,c,scale) end return f end function laweffect(yr,rate,start) if yr >= start 1.0-exp(-rate*(yr-start)) else 0.0 end end statecoefpert = log(1.25)/2.0 include("model.jl") ## load the model, this is a separate file so we can save the output unless it changes modeljoined = @chain leftjoin(joined,stategroupdf; on = :STUSAB, makeunique=true) begin @subset(.! ismissing.(:crudedeathrate)) end lawdates = @chain leftjoin(DataFrame(STATE=1:55),concarry; on=:STATE) begin @subset(:STATE .< 60) @rtransform(:year = ismissing(:year) ? 3000 : :year) @orderby(:STATE) end centers = [i for i in (minimum(modeljoined.year)-4):4:(maximum(modeljoined.year)+4)] width = 7.0 modl = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup), centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),lawdates.year,statecoefpert) setadbackend(:reversediff) Turing.setrdcache(true) savedfile = "./saved/samples.dat" global s = [] if stat(savedfile).mtime > stat("./model.jl").mtime global s = deserialize(savedfile) else global s = sample(modl,NUTS(500,0.8),MCMCThreads(),500,3) serialize(savedfile,s) end
Having sampled the model, we can then compute summary graphs of the results. We will show the posterior density for the coefficient of the magnitude of the law effect for each state separately.
Also we will plot the actual log(homicide rate) together with the state specific counterfactual estimates.
lawco = group(s,"statelawcoef") statenames = Dict() for (st,stusab) in Iterators.zip(lawdates.STATE,lawdates.STUSAB) statenames[st] = stusab end function plotlawts(s,lawdates,modeljoined) let pl = [] lawco = group(s,"statelawcoef") for i in @select(@subset(lawdates,:year .< 2022),:STATE).STATE den = density(s[:,Symbol("meanlawcoef"),:] .+ lawco[:,Symbol("statelawcoef[$i]"),:], title = "State $i = $(statenames[i])", legend=false,xlim=(-log(2),log(2))) plot!([(0.0,0.0),(0.0,5.0)], color="red",ylim=(0.0,5.0)) push!(pl,den) end println("State law coefficient values") display(plot(pl...; size = (1000,1000))) display(density(s[:,:lawrate,:], title= "Law effect onset rate (1/yr)")) modeljoined.logdrate = log.(modeljoined.crudedeathrate) regions = Dict() for r in eachrow(stategroupdf) regions[r.STATE] = (group=r.stgroup,code=r.STUSAB) end pl = [] samps = sample(1:500,10) for st in unique(modeljoined.STATE) sub = @subset(modeljoined, :STATE .== st) region = regions[st].group p = plot(sub.year,sub.logdrate,linewidth=3,title="$(st) = $(regions[st].code)",ylim=(-0.5,2.75)) push!(pl,p) for samp in samps statelawcoef = s[samp,Symbol("statelawcoef[$st]"),1] statecoefs = [s[samp,Symbol("statecoefs[$st][$i]"),1] for i in 1:length(centers)] statebase = s[samp,Symbol("statebase[$st]"),1] regioncoefs = [s[samp,Symbol("regioncoefs[$region][$i]"),1] for i in 1:length(centers)] lawrate = s[samp,:lawrate,1] startlaw = lawdates[lawdates.STATE .== st,:].year[1] startlaw = ismissing(startlaw) ? 3000 : startlaw years = 1979:2020 pred = [statebase + timeser(yr,regioncoefs .+ (statecoefpert .* statecoefs),centers,width) for yr in years] plot!(years,pred; color="orange",alpha=0.5) if startlaw < 3000 plot!([(startlaw,-10.0),(startlaw,10.0)],color="red") end regpred = [statebase + timeser(yr,regioncoefs ,centers,width) for yr in years] end end println("Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)") display(plot(pl...; size =(2000,3000), legend=false)) end density(s[:,:meanlawcoef,:],title="Mean law coefficient") end plotlawts(s,lawdates,modeljoined)
State law coefficient values Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)
One plausible explanation for the relative lack of clear measurable effect from permitless carry laws is that permitless actually may have had very little effect on how many people carry firearms in public. Many of these states were shall-issue and many people who wanted to carry may have had CCW permits issued before the permitless law passed. Therefore the change in number of people carrying may have been minimal.
Using the date on which states converted to shall-issue permits, with a similar analysis, we can determine whether more of an effect was visible after that legal change. Some states converted from no-issue to shall issue, while others converted from may-issue to shall-issue.
# dates collected from https://en.wikipedia.org/wiki/History_of_concealed_carry_in_the_United_States shallissdates = CSV.read("data/ShallIssueDates.csv",DataFrame) shallissue = leftjoin(fipscodes,shallissdates; on = :STUSAB) silawdates = @chain leftjoin(DataFrame(STATE=1:55),shallissue; on=:STATE) begin @subset(:STATE .< 60) @rtransform(:year = ismissing(:ShallIssueDate) ? 3000 : :ShallIssueDate) @orderby(:STATE) end modl2 = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup), centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),silawdates.year,statecoefpert) savedfile = "./saved/shallisssamples.dat" global sisam = [] if stat(savedfile).mtime > stat("./model.jl").mtime global sisam = deserialize(savedfile) else global sisam = sample(modl2,NUTS(500,0.8),MCMCThreads(),500,3) serialize(savedfile,sisam) end
plotlawts(sisam, silawdates,modeljoined)
State law coefficient values Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)
Many states that have passed constitutional carry laws have done so recently enough that there is no data on homicides available from CDC Wonder yet.
For those which passed the law long enough ago to have some post-law data, still many of them were in the last few years leaving only a few years of post-law data. The last few years is a time period where homicides have increased nearly nationwide and it is challenging to estimate the effect of a law when it is imposed on top of a nonlinearly changing trend. Estimates of a law effect created by comparing the actual data to counterfactuals estimated as small perturbations to the regional trend lead to estimates of the average law coefficient that has uncertain sign, with a most likely value of about 0.05 or so, meaning that gun homicides may increase by a factor of around 1.05, however the ability to resolve this number is so poor that the number could be anything from about -0.1 to +0.25 (multiplying the crime rate by anywhere from 0.9 to 1.28.
Any kind of convincing evidence for a consistent average net positive or negative effect of passing constitutional carry laws simply isn't there. What's more, we have decent bounds on the size of this effect for most states, it should be somewhere between perhaps -0.25 and 0.25 on the log scale, meaning approximately the rate of homicides would be multiplied by some number in the range 0.78 to 1.28.
Because many shall-issue laws were passed in the 1990's there is much better evidence for their long term effect on crime. However, using the same methodology as above the average effect is essentially symmetrically distributed around 0.0 with a posterior credible range of about -0.1 to 0.1. The only reasonably strong evidence we have in this data for a real effect is that Texas may have reduced crime with a coefficient around -0.3 and perhaps Pennsylvania caused an increase around 0.2, with Wisconsin possibly around 0.2. The PA estimate does see the counterfactual diverge around the time of the law, but the Wisconsin counterfactual diverges 20 years before the law, suggesting that that effect is really an artifact of a poor counterfactual estimate.
Firearms sales are driven by multiple factors. One is the fear of "bans" which drove a lot of sales around the time of the Sandy Hook school shooting in Dec 2012. And in general sales have increased over the last decade including a sudden increase after the turmoil of 2020, with women and racial minorities making up a historically unprecedented fraction of new sales (close to 50%).
Firearms sales can plausibly be both a cause of firearms crime (when people purchase firearms to commit crime), a response to violence (when people purchase firearms for defense after periods of crime), and even in some cases a deterrent to crime (when criminals meet with greater defensive responses than before). It is hard to say which of these effects dominates, and in fact most likely all of these occur at different times and different places. One thing is clear however, and that is that it is entirely possible for gun ownership to increase while crime decreases, or for crime to increase while gun ownership per capita decreases. Those who somehow wish to draw a direct causal link in the US between more guns and more gun crime simply do not have an easy obvious convincing argument.
The evidence suggests that in the post NYSRPA v Bruen era, forcing CA and NY and others to shall-issue their permits will by itself have very little effect on overall gun crime. More likely if there is an increase it will be due to the turmoil undergoing the whole country at the moment, with high levels of inflation and the war in Ukraine driving food prices up and crime in general increasing over the last few years.